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Abstract 

The coupled cluster method (CCM) is a powerful and widely applied technique of modern-day 
quantum many-body theory. It has been used with great success in order to understand the 
properties of quantum magnets at zero temperature. This is due largely to the application of 
computational techniques that allow the method to be applied to high orders of approximation 
using localised approximation schemes, e.g., such as the LSUBm scheme. In this article, the high- 
order CCM formalism for the ground and excited states of quantum magnetic systems are extended 
to those with spin quantum number s > |. Solution strategies for the ket- and bra-state equations 
are also considered. Aspects of extrapolation of CCM expectation values are discussed and future 
topics regarding extrapolations are presented. 
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I. INTRODUCTION 



The coupled cluster method (CCM) 

iiB.aflfly.flaa 

is a well-known method of 

quantum many-body theory (QMBT). The CCM has been applied with much success over 
the last fifteen or so years in order to s tudy qu antum magnetic systems at zero temperature 



use of computer-algebraic implementations 




52l | ). In particular, the 
32| of the CCM for quantum systems of 



large or infinite numbers of particles has largely been found to be very effective with respect 
to these spin-lattice problems. This approach uses localised approximation schemes, such 
as the LSUBm approximation. For the LSUBm scheme, the extent of the locale over which 
multi-spin correlations are explicitly included in the approximation is defined by the index 
m. The ground- and excited-state expectation values are often extrapolated in the limit 
m -H> oo. In this article we focus on the development of new high-order CCM formalism for 
the ground and excited states of lattice quantum spin systems with spin quantum number 



s > 



The solution the ket- and bra-state equations is also considered. Various aspects 



of the extrapolation of CCM expectation values are considered and future topics regarding; 
extrapolations are described. The high-order CCCM code is freely available online 36 1. 



II. CCM GROUND-STATE FORMALISM 

The ket and bra ground-state energy eigenvectors, and of a general many-body 
system described by a Hamiltonian H, are given by 

H\^) = Eg\^); {^\H = Eg{^\. (1) 

Furthermore, the ket and bra states are parametrised within the single-reference CCM as 
follows: 

(§| = ($|5e-^ ; S = l + J2SiCT. (2) 

It may be proven from Eqs. ([1]) and ([2]) in a straightforward manner that the ket- and 
bra-state equations are thus given by 

($|C7e-^i/e^|$) =0, V/ 7^ ; (3) 
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mSe-^[H,Ct]e^\<!>) =0, ^ . (4) 

The index / refers to a particular choice of cluster from the set of (Np) fundamental clusters 
that are distinct under the symmetries of the crystallographic lattice and the Hamiltonian 
and for a given approximation scheme at a given level of approximation. We note that these 
equations are equivalent to the minimization of the expectation value of H = {^\H\'$) with 
respect to the the CCM bra- and ket-state correlation coefficients {Si, Si}. We note that Eq. 
([3]) is equivalent to 6H/6S1 = 0, whereas Eq. (jlj) is equivalent to 6H/6S1 = 0. Furthermore, 
we note that Eq. ([3]) leads directly to simple form for the ground-state energy given by 

Eg = Eg{{Si}) = (<l>|e-^i/e^|$) . (5) 

The full set {Si, Si} provides a complete description of the ground state. For instance, an 
arbitrary operator A will have a ground-state expectation value given as 

A = {^A\^) = mSe-^Ae^\<^) = A [{Si, Si}) . (6) 

The similarity transform of A is given by, 

A ^ e-'Ae' = A + [A, S] + ^[[A, S],S] + --- . (7) 

Finally, we remark that the CCM provides exact results in the limit of inclusion of all possible 
clusters in 5* and S. However, this problem is often impossible to solve in a practical sense. 
Hence, we generally make approximations in both S and S. The three most commonly 
employed approximation schemes previously utilised have been: (1) the SUBn scheme, in 
which all correlations involving only n or fewer spins are retained, but no further restriction is 
made concerning their spatial separation on the lattice; (2) the SUBn-m sub-approximation, 
in which all SUBn correlations spanning a range of no more than m adjacent lattice sites 
are retained; and (3) the localised LSUBm scheme, in which all multi-spin correlations over 
all distinct locales on the lattice defined by m or fewer contiguous sites are retained. 

III. HIGH-ORDER GROUND-STATE OPERATORS AND COMMUTATIONS 

We begin the treatment of high-order CCM by introducing the ket-state correlation op- 
erator given, as usual, by 

^ = E E ^n,-,H< ■ (8) 



However, it is important point to note that each of the indices {ii,i2, ■ ■ ■ yk} runs over all 
lattice sites. Furthermore, we assume that there are (/!) orderings of these indices (even for 



s > 



2n 



although we never need to work out these factors explicitly in practice. The index 



/ corresponds to one of the choices of {ii, ■ ■ ■ , «/} for the fundamental set of configurations. 
We may now write a set of high-order CCM ket-state operators, given by 



Nkn.np = El>3 l{l-m-2){l-3)S, 



s7. 



(9) 



The indices k, m, n, and p depend on those sums in the Hamiltonian or of another given 
operator. We note that = ± is^, [s^, s^] = ±s^, and = —2s^. Hence, the 

following commutation relations may be proven: 



slS] 


= Fkst , 




= ~'2Fksl — GkkSk 1 




— r 

^km^k ' 








'^FfnG fi^jjiS 1^ , 




= ~'2GkmS\ — MkkmSt ; 


Sk.Fl] 


= -'^Gl^sl - 2FmMkkmSk - 




-^kmnpS k ; 




= IMkmn^k ■!^kkmnS~k ■ 



'fc ) 



(10) 



We may now write the similarity-transformed expressions of the single-spin operators 
{+,-,2}, as 



s ; a 



Si e" 



e ^sle^ = si 



Si. 



si + Fkst 

S-, - 2Fksl - Gkkst - Fisi 



(11) 



We see that there is a repeated index in Gkk in the similarity transformed version of s . 
Clearly, this term contributes only for systems with spin quantum number s > ^. 
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IV. DERIVING AND SOLVING THE CCM GROUND-STATE EQUATIONS 



We now wish to determine and solve the CCM ket-state equations, where the I-th such 
equation is given by 

Ej = {(^\CYe-^He^\<^) = ,V/^0 . (12) 

(Note that we assume that ($|C7C/|$) = 1 in the above equation). Specific terms in the 
Hamiltonian are now exphcitly written in terms of the high-order CCM operators as: 



TERM 1 


: s'r- 


t 1 Jit ''11 ''J 1 7 J 1 1 


TF,RM 9 




— bj -t- -Tji ■ bj 


TERM 3 




= -2F,sls] - 2G.,sts] - G,,spi - M,,,sts+ - 2F,F,sts] 


TERM 4 






TERM 5 


■.s-s] 


= ~2Fisls] - 2GijSpl - Gus+s] - Mujsts+ - 2F,FjSpl 

^J^i^iyb^ bj rjKjr^ib^ bj -^j-Tj b^ bj -Tj b^ bj 


TERM 6 : 




— — 2FjSj Sj — GjjS^ Sj — Fj Sj Sj 


TERM 7 : 




— z,i '-Jiibi bj 1 j bj 


TERM 8 : 






TERM 9 : 


^7^7 


= AG^jsls] + 2Miijsts] + 2Mij,s+sl + iV«i,s+s+ 



+2G2.S+S+ + 2F,M,ijsts+ + 4F,-G'i,-s+s^ + AF.Fjsls] 
+AF,Gi,sts] + 2F,G,-,4s,^ + 2FiM,nsts+ + 4F,FjG,,sts+ 
+2F,Ffspt + 2FjGusts] + + Fpusfs^ 



+2F-F,.+s^ + F-G,,.+4 + F-F/.+4 
TERM 10 : = + FiS+ 
TERM 11 : ~sj = -2F,s,^ - G^sj - {F^fsf 
TERM 12 : (.7)^ = {slf + 2F..+.,^ + G^f + Hsf? + i^^(^.+)^ 

TERM 13 : s+ = s+ (13) 

(Note that = is implicitly assumed in Eq. fll3p above.) We now "pattern-match" 

the G~ operators to those the relevant terms in the Hamiltonian from Eq. (fT3l) above in 
order to form the CCM equations F/ = of Eq. (fT2|) at a given level of approximation. 



We now define the following new set of CCM bra-state correlation coefficients given by 
xi = Si and x/ = Nb/N{1\)viSi and we assume again that {(^\CjCf\^) = 1. Note that 
Nb is the number of Bravais lattice sites. Note also that for a given cluster I then i/j 
is a symmetry factor which is dependent purely on the point-group symmetries (and not 
the translational symmetries) of the crystallographic lattice and that / is the number of 
spin operators. We note that the factors uj, N, Nb, and (/!) never need to be explicitly 
determined. The CCM bra-state operator may thus be rewritten as 

Nf 

S = l + Nj2S:iCT , (14) 

7=1 

such that we have a particularly simple form for H, given by 

Nf 

H = Nj2S:jEj , (15) 
/=o 

where xq = 1. We note that the Eq is defined by Eq = j^{^\e^^ He^\^) (and, thus, Eq = 
j^Eg) and that Ei is the J-th CCM ket-state equation defined by Eq. ([12]). The CCM ket- 
state equations are easily re-derived by taking the partial derivative oi H /N with respect to 
x/, where 

We now take the partial derivative of H /N with respect to xj such that the bra-state 
equations take on a particularly simple form, given by 

Nf 



A.o)^^, . (17) 

OXj OXj JZi OXi 

The coupled non-linear equations for the ket state Ef = are solved readily, e.g, by using the 
Newton-Raphson method, in order to find the coefficients {xj}. By contrast, the equations 
for the bra state Ej = are easily solved via LU decomposition, although this may only 
be carried out once the CCM ket-state equations have been determined and solved. The 
numerical values of the coefficients {xj} may thus be obtained. We note that this approach 
greatly simplifies the task of determining the bra-state equations because we infer the bra- 
state equations directly from those of the ket-state equations via Eq. (1171) . Thus, we never 
need to evaluate Eq. (jlj) exphcitly. 

We may also solve the ket- and bra-state equations (i.e., Ej = and Ej = 0, respectively) 
via direct iteration. For the case of the ket-state equations this is slightly more complicated 
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because there are non-linear terms with respect the ket-state correlation coefficients {x/}. 
We rearrange the ket-state equations such that the linear terms for xj for the z*^ CCM ket- 
state equation are on the left of the new equation and all other terms are on the right. The 
right-hand side of this new equation is denoted by Ej after dividing through by the factor 
on the left-hand-side for the ket-state correlation coefficients. We may carry out exactly the 
same procedure for the bra-state in order to find E'j, although the problem is linear with 
respect to {xj} in this case. These equations are thus rewritten conveniently for the ket 
state as 

Xj = X2, • • • , Xj-i, • • ■ , XaTj^, x^, ■ ■ ■ , x^, ■ • ■ , x^^) , (18) 

and for the bra state as, 

xj = Ej(xi, X2, ■ ■ ■ , Xj-i, Xj^i, ■ ■ ■ , xnp ; xi, X2, ■ ■ ■ , xn^, xf, X2, ■ ■ ■ , xf, ■ ■ ■ , x\^) .(19) 

Clearly, these equations may be solved for xj and xj by iterating them "directly" until 
convergence. Indeed, the local memory usage is vastly reduced because we do not need to 
store any Jacobian or other large matrix that scales in size with Np. This simple "brute 
force" approach of direct iteration has actually been found to be surprisingly successful. 
However, in practice, it needs to be implemented for large numbers of CPUs used in parallel 
for very large numbers of clusters (e.g., 10^) used in S and S. Clearly, more sophisticated 
solvers for the bra- and ket-state equations that do not demand the memory requirements of 
Newton-Raphson for the ket state and LU decomposition for the bra state and are quicker 
than direct iteration may be implemented. However, this remains a task for the future. 

V. THE EXCITED-STATE FORMALISM 

We now consider how the excited state may be treated using the CCM via a high-order 
approach. We begin by remarking that the excited-state wave function is given by 

l^e) = e^|$) . (20) 

The Schrodinger equation, -Eel^e) = -f^l^e) and the equivalent equation for the ground state 
lead (after some simple algebra) to 

eeX'^l^) =e-^[H,X']e^\^) {=R\^)) , (21) 
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where ee = — Eg is the excitation energy. We note that the excited-state correlation 
operator is written as, 

X^ = Y,'^iCt , (22) 
Equation ( l22ll imphes the overlap relation 

($|^e) = . (23) 

We may now form the basic equations for the excited state, given by 

e^X^ = {(!>\CYe'^[H,X^]e^\^) ^1^0 , (24) 

which is a generalized set of eigenvalue equations with eigenvalues eg and corresponding 
eigenvectors Xj. We note that the choice of clusters for the excited-state may be different 
from those for the ground state. For example, the ground state for the Heisenberg model 
on bipartite lattices is in the subspace = = O5 whereas the excited state has 

= = +1- The number of excited-state "fundamental" clusters that are distinct 

under the translational and point-group symmetries of the lattice and Hamiltonian is given 
by%. 



VI. HIGH-ORDER EXCITED-STATE OPERATORS AND COMMUTATIONS 

In a similar manner as for the ground-state, we now define excited state operator via 

^^ = E E K-,n<---st (25) 

I n, 

where the indices {ii, - ■ ■ , ii} run over all lattice sites. We assume explicitly again that there 
are (/!) orderings of the indices (even for s > |). The index / corresponds to one of the 
choices of {ii, ■ ■ ■ ,ii} for the fundamental set of configurations for the excited state, such 
that Xf = ... j^. We now also define the further high-order operators for the excited state, 
given by 

Rkmn = J2l>2 ~ 1)(^ ^ '^k,m,n,i4,--- ,ii ' ' ' 

Tkmnp = J2l>3 J2i5,-,ii ~ 1)(^ ~ 2)(/ — 'i)Xk,m,n,p,i5,--,ii ^ts ' ' ' 
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The following commutation relations may also be proven: 



[si P. 



i^k 5 Qmn] 



Pkst , 

~'^PkSl — QkkSk , 
Qkm^k 1 
Pkmn^k 1 
'^PmQkm'^k ' 

'^Qkm^k Pkkm^k ' 

'^Q\mSi '^PmRkkm^k ^^PrnQkm^k 
Tkmnp'^k ' 



(27) 



VII. DERIVING AND SOLVING THE EXCITED STATE EQUATIONS 



We now wish to determine and solve the CCM excited-state equations given by Eq. ([24 
Specific terms in the Hamiltonian are now explicitly written in terms of the new excited-state 
high-order CCM operators as: 



TERM 1 : e-^[sls],X^]e^ 
TERM 2 : e-^[sls^ , X^]e^ 
TERM 3 : e-^[s,"s-, X"]e^ 



TERM 4 : e-^[s+s", X"]e 



TERM 5 : e-^[s-s", X"]e^ 



TERM 6 : e-^fs+s", X'=]e^ 
TERM 7 : e'^lsi , X^]e^ 
TERM 8 : e'^lst , X^-^ 
TERM 9 : e-^[s- , X'']e^ 



Pisfs'^ + P^F^sfs^ + Pjspl + PjFists+ + QijS+s] 



-2P,F,sts] - P.G,,sts^ - P^Ffsts| - 2P,sls] - 2P,F,sp 

2PjFiS,^ Sj 2PjGijS^ Sj 2PjF{FjSj^ Sj 2C^j^jS^ Sj 
—2QijFjS^ Sj — QjjSj Sj — QjjFiSj^ Sj — RijjSi Sj 

-2PjF,s+sl - PjGiists+ - PjFfspt - 2P^sls] - 2PiFjSp 
-2P,F,sts] - 2P,G,,stsj - 2PiFiFjsts^ - 2Q,jSpl 

2QijFiSj^ Sj QiiS^ Sj QiiFjS^ Sj RiijS^ Sj 
—2Pjsfsj — Qjjsf s'j — 2PjFjsf s'j 

2PiSj Sj QiiSj^ Sj 2fjPjSj Sj 



APiFjsls] + 4PjGj,s+s^" + 2PiGjjSpl + 2PiMijjsts^ 
+APiFiF,sts] + APiFjGi.sts^ + 2P,FjG,jS+s+ + 2P,FiF]st 
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+AQijs'iS'j + AQijFjSp'i + AQijFisfs] + AQijdjsts^ 
+4Q,,FiFjsts+ + 2Riijsts] + 2Ri,jFjsts+ + 2Rij,spt 

+AP,Fisls'j + APjFiF.spt + '^PjGiists] + 2P,GuF,sts^ 
+2P,FMs] + 2P,F^F,stsl + 4P,F,G,,stsj + 2Q,,F,spi 
+QnG.^4st + QnF^sts^ 
TERM 10 : e-^[s^,X"]e^ = Pisf 

TERM 11 : e-^K-,X^]e^ = -2P,sl - Qusf - 2P,F,st 
TERM l2:e-'[{sff,X']e' = 2P,st si + Q,,{stf + P,{stf + 2P,F,{stf 

TERM 13 : e-^[s+, X'^je^ = (28) 

(Note that s~|$) = is again implicitly assumed in Eq. 0281) above.) Again, we now 
"pattern-match" the Gj operators (this time with respect to the fundamental set of the 
clusters in the excited state) to those the relevant terms in the Hamiltonian from Eq. (!24|l 
above in order to form the CCM excited-state equations at a given level of approximation. 
By contrast to the case for the ground state, we see that the high-order operators of Eq. 
fl2^ are in linear in those terms in Eq. fl2S]) . We choose the eigenvalue of lowest value to be 
our result, and this method was found to provide good results in regions of the parameter 
space for which the model state was a good choice. Again we note that we have formed 
an eigenvalue problem, which is readily solved using a standard eigenvalue solver. However, 
the computational problem thus formed uses local memory that scales with the number of 
fundamental clusters used in the excited state, i.e., as Nj^. 

The eigenvalue equations of Eq. (12^ may be iterated directly in order to solve them. 
We denote the matrix for the eigenvalue problem of Eq. ([21]) by B and we denote the 
eigenvectors hy y = (Aff , ■ ■ ■ , ^ . Hence, we iterate directly the eigenvalue equation 
given by 

By = Xy . (29) 

This is just the well-known "power iteration" method and the ratios of Xj in successive 
iterations yields the relevant eigenvalue. However, the eigenvalue determined in this manner 
is the eigenvalue of largest magnitude, Amax, rather than the lowest (generally the one of 
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smallest magnitude Amin for our purposes) that we wish to obtain here. Thus, we use find 
the eigenvalue of smallest magnitude by using the "shifted" power iteration method. Once 
Amax has been found, we then solve the following eigenvalue equation by direct iteration: 

(5 - AMAx/)y' = AV . (30) 

This process ought to converge to an eigenvalue A' = Amin — Amax- Indeed, this was found to 
be the case for those spin models for which the model state was a "good choice" . Indeed, the 
lowest-valued eigenvalue obtained in this manner agreed perfectly with those results for the 
eigenvalue of lowest values obtained via a complete diagonalization of the matrix eigenvalue 
problem at every level of approximation. 



VIII. EXTRAPOLATION OF EXPECTATION VALUES 

In practice, we often need to extrapolate individual LSUBm or SUBm-m (etc.) expec- 
tation values y{m) in the limit m — oo. Indeed, the extrapolation of COM expectation 
values is pivotal to its practical use for quantum magnetic systems. However, there are 
many aspects of such extrapolations that we still do not understand and that have not been 
tested. Furthermore, the extrapolation of CCM expectation values poses two distinct prob- 
lems: a) no provable exact rules of extrapolation are known (unlike, e.g., finite-size exact 
diagonalisations); and, b) we have only small numbers of data points (i.e., 7 or 8 at most for 
ID systems). Hence, until now, only heuristic or "ad hoc" schemes have been used. These 
heuristic schemes are "parametric" in the sense that the data is fitted to specific functions 
(e.g., polynomials) that are defined prior to fitting the data and where the coefficients of 
these functions are the basic "parameters." However, one should note that the results of 
such parametric/ad hoc extrapolation procedures have been shown to yield consistently valid 



results for an extremely wide range of quantum spin systems, see, e.g., Refs. [22|, |26|, l32 |. 
Indeed, the use of such schemes has been shown to provide improved results in all cases 
where the model state might even be suspected to be a reasonable starting point. This is a 
very strong vindication of the use of these simple extrapolation procedures. 

We now turn to specific instances of such extrapolation schemes, and we start by noting 
that a common scheme for extrapolating the ground-state energy is given by 

2/("^)scHEME 1 = ao + aim~^ + b2m~^ . (31) 
11 



Thus far, least-squares fits of the data to this scaling rule have been performed, where the 
extrapolated value y{m — > cxd)scheme i is given by ai. For other expectation values, such 
as the sublattice magnetisation and excitation energy gap, a common scheme is given by 

yini)scuEME 2 = bo + bim~^ + bi'm''^ . (32) 

In this case, the extrapolated value y{m oo)scheme 2 is given by bo. Finally, another 
common scheme is given by 

y{m)scuEME 3 = Co + cim"^ . (33) 

The extrapolated value y{m — > cxd)scheme 3 is given by cq. Other schemes have included 
Fade approximant s [26 1 and similar schemes to Eqs. 0311) and (1321) for fixed, though non- 



integer, exponents 



481 ] have been used. For example, the sublattice magnetisation 
in Refs. 45|, |46|, |47|, |48| was found to scale as: y{m) = cIq + dim^^'^ + d2m"^'^. 



In future, we would ideally like to establish exact rules of scaling of expectation values 
with approximation level and/or go to much higher orders of approximation so that we have 
more data points to extrapolate. Evidence for exact rules of scaling is supported by the fact 
that the ground-state energy of the Heisenberg model on bipartite lattices such as the linear 
chain and square lattice appear to follow the scheme of Eq. ( l3Ti) quite well. By contrast, 
their magnetisations and excitation energy gaps appear to follow the scheme of Eq. (!32l) . 
Few other general rules seem to exist and certainly no such rules have, as yet, been proven 
mathematically. Furthermore, we are must use computationally intensive methods in parallel 
to solve for levels of approximation currently available. Further increases in approximation 
level might be possible in future due to computational improvements, although we might still 
be restricted to fairly small numbers of data points for the current approximation schemes. 
Hence, both of these goals may be not be achieveable in practice. In the absence of either of 
these goals, one might wish to extrapolate using a variety of extrapolation schemes in order 
to determine (in broad terms only) the amount of their mutual agreement. Furthermore, 
odd and even series of (e.g., LSUBm expectation values) ought to extrapolate to the same 



value (see Ref. 50j). Again, this might yield a rough idea of errors of extrapolation. 

Future research might also concentrate on the establishment of new approximation 
schemes that do not scale as quickly as LSUBm and SUBm-m schemes and yet still appear 
scalable with approximation level. However, it is likely that the number of fundamental 
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clusters in any such scheme will still increase exponentially with approximation level. It 
would be beneficial to carry out a fully "statistical" exploration of the topic of parametric 
extrapolations such as those of Eqs. (1311) to (1331) based on only small numbers of data 
points. This might yield better extrapolations, but perhaps more importantly it would hope 
to establish with more mathematical rigour the estimate of the error in the extrapolated 
values. Until now, the extrapolated figures have generally been presented with no estimated 
degree of error. Furthermore, this "statistical" approach might also include an exploration 
of, e.g., least-squares, weighted least-square and maximum likelihood methods fitting meth- 
ods of the parametric scaling laws to the data and the effect that this would have on the 
extrapolated results. Indeed, hitherto, only ad hoc/parametric models or rules of the types 
shown in Eqs. (13 ip to ([3SD have been used to carry out extrapolation of CCM data using 
least-squared methods. Finally, a outline of "best practice" for extrapolating CCM data 
(i.e., which extrapolation schemes and methods to use and exactly what to report) for given 
types of expectation values might prove very useful. These topics remain for future study 
however. 
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